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, Abstract 

In this paper we formulate a geometric theory of the mechanics of growing solids. Bulk growth is modeled 
q ■ by a material manifold with an evolving metric. Time dependence of metric represents the evolution of the 

y-^- ' stress-free (natural) configuration of the body in response to changes in mass density and "shape" . We show 

that time dependency of material metric will affect the energy balance and the entropy production inequality; 
both the energy balance and the entropy production inequality have to be modified. We then obtain the 
governing equations covariantly by postulating invariance of energy balance under time-dependent spatial 
diffeomorphisms. We use the principle of maximum entropy production in deriving an evolution equation for 

Othe material metric. In the case of isotropic growth, we find those growth distributions that do not result in 
residual stresses. We then look at Lagrangian field theory of growing elastic solids. We will use the Lagrange- 
|_} ■ d'Alembert's principle with Rayleigh's dissipation functions to derive all the governing equations. We make 

y—n ' an explicit connection between our geometric theory and the conventional multiplicative decomposition of 

w ■ deformation gradient F = F e F 9 into growth and elastic parts. We linearize the nonlinear theory and derive a 

linearized theory of growth mechanics. Finally, we obtain the stress-free growth distributions in the linearized 
M | theory. 
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1 Introduction 

In classical continuum mechanics, one usually models mass-conserving bodies. The traditional framework of 
continuum mechanics is suitable for many practical applications. However, in some natural phenomena mass is 
added or lost in a deformation process. This is particularly important in biological systems where growth and 
remodeling are closely linked to mechanical loads. In the case of soft tissues, elastic deformations are large and 
theory of linear elasticity is not adequate. This has been the main motivation for the extensive studies of large 
deformations in biomechanics in recent years (see Cowin and Hegedus [1976]; Skalak, et al. [1982]; Fung [1983]; 
Naumov [1994]; Hoger [1997]; Humphrey [2003]; Klarbring, et al. [2007] and references therein). 

It has been realized for a long time that mechanical forces directly affect growth and remodeling in biological 
systems [Hsu, 1968]. A continuum theory capable of modeling biological tissues must take into account changes 
of mass and the coupling between growth/remodeling and mechanical stresses. In continuum mechanics, one 
starts by postulating that a body is made of a large number of "material points" that can be treated as 
mathematical points. A material point consists of "enough" number of particles (atoms, molecules, cells, etc.) 
such that it can represent the mechanical properties of the body, e.g. density. Material points are then identified 
with their positions in the so-called reference configuration. This is called the material manifold. It is always 
assumed that the body is macroscopically stress-free in the material manifold. The material manifold is not 
necessarily Euclidean and even not Ricmannian; in general, material manifold is a Riemann-Cartan manifold in 
the case of solid bodies with distributed dislocations, for example. It is relevant to mention that in most of the 
existing formulations of finite-strain plasticity, instead of working with a Riemann-Cartan material manifold, one 
assumes a multiplicative decomposition of the deformation gradient into clastic and plastic parts, i.e. [Eckart, 
1948; Kroner, 1960; Lee and Liu, 1967; Lee, 1969] 

F = Tip = F e F p , (1.1) 

where ip is the deformation mapping. This means that locally the material deforms plastically followed by elastic 
deformations to ensure compatibility. In other words, one assumes that both the material and the ambient space 
manifolds are Euclidean and one locally decomposes the total deformation mapping into incompatible elastic and 
plastic parts. Motivated by plasticity, in the case of growth several researchers [Kondaurov and Nikitin, 1987; 
Takamizawa and Matsuda, 1990; Takamizawa, 1991; Rodriguez, et al., 1994] postulated a similar decomposition 
of F into elastic and growth parts, i.e. 

F = F e F s . (1.2) 

This has been the fundamental kincmatical idea of all the existing models of growth mechanics to this date (see 
Ben Amar and Goriely [2005]; Hoger [1997]; Lubarda and Hoger [2002] and references therein). 

Recently, Ozakin and Yavari [2009] introduced a geometric theory of thermoelasticity in which thermal 
strains are buried in a temperature-dependent Ricmannian material manifold. In that theory a change of 
temperature leads to a rescaling of the material metric with a clear physical meaning. In this paper we introduce 
a geometric theory for growing bodies using similar ideas. We should mention that the analogy between growth 
and thermal distortions was first realized by Skalak and his coworkers [Skalak, et al., 1996]. 

There are two possibilities in a growth process: (i) the number of material points is fixed, and (ii) material 
points are removed or are added. Note that in a continuum model material points are assumed to contain 
several (or a large number of) particles, cells, etc. Erosion or accretion of cells corresponds to changes in volume 
(and the corresponding mass) and shape of the material body. In our continuum model, similar to many of 
the earlier models on bulk growth, we assume that the number of material points is fixed. This means that we 
work with a fixed set B as the material manifold and model growth by allowing B to have an evolving geometry. 
Consider a two-dimensional problem, where the relaxed state of the material is described by a surface. If the 
bulk of the material grows as, for example, in a thin shell of biological material undergoing cell division, the 
shape of the surface describing the relaxed state will change. The stresses for a given configuration should be 
calculated in terms of the map from the surface describing the relaxed state, to the current configuration. 

While the multiplicative decomposition of the deformation gradient has been a source of useful approaches 
to nonlinear problems, we believe that in many cases, such an approach obfuscates the underlying natural 
geometry. A multiplicative decomposition seems natural if one starts with a stress-free material body and 
considers processes such as plasticity which, in general, induce stresses. However, an initial stress-free Euclidean 
configuration may not even exist in certain problems. Mathematically, one can still consider an incompatible 
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local deformation that brings the material to a relaxed, Euclidean state, and measure deformations from this 
state, as in the multiplicative decomposition described above. However, we believe a more natural way of 
looking at the problem involves treating the material manifold as a non-Euclidean manifold, and giving its 
geometry explicitly in terms of the physics of this problem. In passing we should mention Miehe's [Miehe, 
1998] work in which instead of F = F e F p , he introduces a "plastic metric", although with no clear physical 
meaning/interpretation for this metric. We will come back to a geometric interpretation of F = F e F 9 in §3. 

In this paper, we model growth by introducing a Riemannian material manifold with an evolving metric. 
As will be seen, this formalism is very similar to the approach of [Ozakin and Yavari, 2009] to thermoelasticity, 
however, there are a few important modifications. First, although we had a version of mass conservation in 
thermoelasticity in terms of the changing material manifold, for the case of growth, mass will in general be 
added to (or removed from) the material body; we will have a mass balance. Thus, one can represent the 
amount of mass being added (or removed) in terms of the changes of the differential form describing material 
mass-density. Secondly, for the case of thermal stresses, the material metric was explicitly given in terms of the 
temperature, but no such simple dependence exists for the material metric in biological growth. We will begin 
by exploring the consequences of various simple modes of growth, such as a cylindrically symmetric growth 
represented by a radius-dependent conformal scaling of the metric. Assuming simple constitutive relations, 
we will write the equations for equilibrium configurations in terms of the time-dependent metric, much like 
the case in thermoelasticity. We will also establish the connection to the formulations involving multiplicative 
decomposition of the deformation gradient. 

One should note that there is in fact no guarantee that a time-dependent Riemannian metric and its Levi- 
Civita connection is capable of modeling all kinds of growth. We believe that at the very least one needs to 
consider time-dependent connections with torsion, however, "Riemannian growth" is a good starting point. We 
aim to investigate the case of growth with torsion, as well as non-metricity in future communications. 

Efrati, et al. [2009] have recently studied similar problems in the framework of linearized elasticity by mod- 
ifying the definition of linearized strain. Here, we start with nonlinear elasticity and instead of modifying any 
definition of strain will work with an evolving material manifold. We should mention that the idea of using 
differential geometry in elasticity goes back to more than fifty years ago in the work of Eckart [1948] who realized 
that the stress-free configuration of a material body evolves in time and an Euclidean stress-free configuration 
is not always possible. Later developments are due to Kondo [Kondo, 1955a,b] and Bilby [Bilby, et al., 1957, 
1955]. 

There have been growth models in the literature using mixture theories. For growth mechanics purposes, a 
mixture theory is certainly more realistic than a mono-phasic continuum theory. However, in this paper for the 
sake of simplicity and clarity of presentation, we restrict ourselves to mono-phasic continua. We should mention 
that our ideas are similar, in spirit, to those of Rajagopal and Srinivasa [2004b] who have been advocating the 
idea of material bodies with evolving natural configurations. Here, we work in a fully geometric framework and 
model a growing body by a continuum that has an evolving Riemannian material metric. One should note that 
this is a very special case of a possible evolving material manifold that we believe is sufficient for bulk growth 
purposes. In particular, we work with Levi-Civita connections that are torsion-free. We should also emphasize 
that we are not, by any means, questioning the usefulness of the traditionally used F = F e F g decomposition 
of deformation gradient. However, we believe that although the existing models based on this multiplicative 
decomposition have been very useful in growth mechanics (see Ben Amar and Goriely [2005] ; Garikipati, et al. 
[2004] for some concrete examples.) they all lack a rigorous mathematical foundation. 

This paper is organized as follows. In §2 we modify the existing geometric theory of elasticity for growing 
solids. We show that energy balance has to be modified and then study its covariance. We study the entropy 
production inequality and the restrictions it imposes on constitutive equations. We also show how Principle 
of Maximum Entropy Production can be used to obtain thermodynamically-consistcnt evolution equations for 
the material metric. We then look at isotropic growth and model it by a time-dependent rcscaling of an initial 
material metric. We solve three examples of isotropic and non-isotropic growth analytically. We then discuss how 
an evolving material manifold can be visualized using embeddings. In the last part of this section we obtain 
stress-free isotropic growth distributions. We then study growth in the Lagrangian field theory of elasticity 
and show how all the governing equations can be obtained using the Lagrangc-d'Alcmbert principle and using 
Rayleigh's dissipation functions. In §3, we make a connection between the exiting theories of growth based on 
the decomposition F = F e F g and the present geometric theory. The nonlinear geometric theory is linearized 
in §4. In particular, we obtain those isotropic growth distributions that are stress free in the linearized setting. 
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Conclusions are given in §5. 



2 Evolving Material Metrics and Bulk Growth 

There have been previous works on continuum mechanics formulation of bodies with variable mass (see Lubarda and Hoger 
[2002]; Ben Amar and Goriely [2005]; DiCarlo and Quiligotti [2002]; Epstein and Maugin [2000], and references 
therein). In these works it is assumed that the growth part of deformation gradient is an unknown tensor field 
and its evolution is given by a kinetic equation. In writing energy balance, the corresponding thermodynamic 
forces show up. In the present geometric theory, we work with an evolving material manifold instead of intro- 
ducing new fields other than material mass density. Before going into the details of the proposed theory, let us 
first briefly review the geometric theory of classical elasticity. 



Geometric Elasticity. A body B is identified with a Ricmannian manifold B and a configuration of B is a 
mapping tp : B — > S, where S is another Ricmannian manifold. The set of all configurations of B is denoted 
by C. A motion is a curve c : R — > C; t i— > tpt in C. It is assumed that the body is stress free in the material 
manifold. For a fixed t, </?t(X) = cp(X,t) and for a fixed X, <^x(£) = </>(X, t), where X is position of material 
points in the undeformed configuration B. The material velocity is the map V t :B->R 3 given by 

V,(X ) =V(X,t) = ^ = ^ x ( t ). (2.1) 
The material acceleration is defined by 

A t (X) = A(X,i) = = |vx(t). (2.2) 

In components 

dV a h 

A a = —+^ c V b V c , (2.3) 

where j£ c is the Christoffcl symbol of the local coordinate chart {x a }. Note that A does not depend on the 
connection coefficients of the material manifold, tpt is assumed to be invcrtible and regular. The spatial velocity 
of a regular motion <p t is defined as 

v t : Vt (B)^M. 3 , v t = V t o VJt - 1 , (2.4) 
and the spatial acceleration is defined as 

a = v = -^ + V v v. (2.5) 
at 

In components 

a =^r + w v + ^ vv - (2 - 6) 

Let tp : B — > S be a C 1 configuration of B in S, where B and S are manifolds. Deformation gradient is the 
tangent map of tp and is denoted by F = Tp. Thus, at each point X € B, it is a linear map 

F(X) : T X B -> T v(x) 5. (2.7) 

If {x a } and {A A } are local coordinate charts on S and B, respectively, the components of F are 

F a A (X) = J|^(X). (2.8) 

Suppose B and S are Riemannian manifolds with inner products ((, )) G and ((, )) based at X € B and xeS, 
respectively. Transpose of F is defined by 

F T : T X S -> T X S, ((FV, v» g = ((V, F T v)) G We T x 6, v e T X S. (2.9) 
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In components 

(F T (X)) A a = g ab ( X )F b B (X)G AB (X). 
The right Cauchy-Green deformation tensor is defined by 

C(X) : T X B -» T*B, C(X) = F(X) T F(X), 

where g and G are metric tensors on S and B, respectively. In components 

C A B = (F T ) A a F a B . 



One can show that 



C b = ^*(g)=F*gF, i.e. C A1 



{g ab o<p)F a A F b B . 



(2.10) 
(2.11) 

(2.12) 
(2.13) 



For bulk growth, we assume that the material manifold B remains unchanged but the metric evolves, i.e. 
G = G(X, t) 1 . When mass is added or removed, the stress-free state of the body changes. Local changes in 
mass change the stress-free configuration of the body. This is modeled by a timc-dcpcndcnt material metric that 
represents local changes in volume and "shape" in the relaxed configuration (see Fig. 2.1). In §3, we will make 
a connection between this approach and the conventional F = F e F p decomposition of deformation gradient. 




(B,G(t)) 



Figure 2.1: Deformation of a growing body as a map between a Riemannian material manifold with a dynamic metric and an 
ambient space with a fixed background metric. 



Incompressibility. In growth mechanics it is usually assumed that elastic deformations are incompressible. 
In the classical theory in which F = F e F 9 is assumed, incompressibility implies detF e = 1, i.e. all the volume 
changes are due to growth. In the geometric theory the following relation holds between volume elements of 
(B,G) and (<S,g): 

dv = J dV, (2.14) 

where 



det 2 

J =Vd^! dctF - ^ 

Incompressibility of elastic deformations means that J = 1. Note that even when J = 1, still dv is time 
dependent as a result of the time evolution of the material metric that makes dV time dependent. In other 
words, an observer in the ambient space sees changes in volume that are only due to volume changes in the 
material manifold. We will show the equivalence of J e = 1 in the classical theory with J = 1 in the geometric 
theory in both some simple examples in §2.9 and in the general case in §3. 



1 In mathematics, evolving metrics have been studied extensively. The most celebrated example is Ricci flow [Hamilton, 1982; 
Topping, 2006] that was used in proving Poincarc Conjecture by Perelman [2002], Interestingly, for a seemingly very different 
application, i.e. growth mechanics, an evolving geometry plays a key role. 
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2.1 Energy Balance. 

Let us look at energy balance for a growing body. The standard material balance of energy for a subset U C B 
reads [Yavari, et al., 2006] 

jti p °{ E+ \ ((v ' v>> ) dv =L p ° (((b ' v>> + r) dv + S m (((t ' v>> + h) dA (2 - i6) 

where E = E(K, N,G,F, g o ip) is the material internal energy density. N, po, B, T, R, and H are specific 
entropy, material mass density, body force per unit undcformcd mass, traction vector, heat supply, and heat 
flux, respectively. 

Wc first note that energy balance should be modified in the case of growing bodies with time-dependent ma- 
terial metrics. Note that when metric is time dependent, material density mass form m(X, t) = /?o(X, t)dV(X., t) 
is time dependent even if p$ is not time dependent. For a subbody IA C B, rate of change of mass reads 



j f Po (K,t)dV(K,t)= [ 



d Po 1 fdG 



dV. (2.17) 



Note that if po is time independent, then the term ^Po^i^f) represents the change in mass due to growth. 
Here, we assume that a scalar field of mass source/sink S m (X.,t) is given. 2 This mass source will change the 
stress- free configuration of the body and (B, G(X, t)) represents the stress-free configuration of the body. 

The rate of change of material metric is a kinematical variable that contributes to power. Therefore, energy 
balance for a growing body with a time-dependent material metric is written as 3 



* Su P °( E+1 2 (<V ' V>> ) ^ = L{ P0 m Y))+R) + P ^--^ 




(((T,V))+H)dA. (2.18) 

dU 



2.2 Covariance of Energy Balance. 



It turns out that in continuum mechanics (and even discrete systems) one can obtain all the balance laws 
using energy balance and postulating its invariance under some groups of transformations. This was intro- 
duced by Green and Rivlin [1964] in the case of Euclidean ambient spaces and was extended to manifolds by 
Marsden and Hughes [1983]. See also Simo and Marsden [1984]; Yavari, et al. [2006]; Yavari and Ozakin [2008]; 
Yavari [2008]; Yavari and Marsden [2009a, b] for applications of covariance ideas in different continuous and 
discrete systems. 

In order to covariantly obtain all the balance laws, we postulate that energy balance is form invariant under 
an arbitrary time-dependent spatial diffcomorphism £ t : S — * S, i.e. 

I L p '° ( e + 5 <<v ' V ' >J ) dV = L {"'" m - v>> + + : tit 



S'„ (e' + i «V, V'))) \dV + J «{T', V')) + H') dA. (2.19) 



2 Note that by definition 



f P0 (X,t)dV= f S m (X,t)dV= f Sm (X,t)dV, 
Ju Ju Ju 



where Sm {X,t) is mass source in the initial material manifold with volume element dV ■ Note also that physically S is given. 

3 Notc that in Lubarda and Hoger [2002] the term analogous to + ^po ^y^fj i s denoted by r g . There, instead of the term 

P°H : they introduce a term pR. g r g . One should note that even if mass is conserved at a point, still a change in shape 
can contribute to energy balance and is captured in our formulation. See also Epstein and Maugin [2000] and Lubarda and Hoger 
[2002]. 
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Note that [Yavari, et al, 2006] 

R' = R, H' = H, p' = Po , T' = &*T, V'=^V + W, 
where W = J^£ t o tp. Note also that 
dG' dG 



G' = G, 



Thus, at t = tc 



dt dt 
d E* 



and E'(X, N', G, F', go^) = E(X, N, G, F, £; g o <p). 



dE <M dE dG 



dE 



£ w g o tp. 



(2.20) 



(2.21) 



(2.22) 



dt dN dt dG dt dgoip 

We also assume that body forces are transformed such that [Marsden and Hughes, 1983] B' A' = £t*(B — A 
Therefore, (2.19) at t = to reads 



dp 1 (dG 



E+\{{V + W,V + W)) 



jf |p (((B, V + W)) + i?) + Po |§ : ^ + 5 m + \ «V + W, V + W))) } dV 
+ [ ({{T,V + W))+H)dA. 

JdU 



(2.23) 



Subtracting (2.18) from (2.23), one obtains 



dpo , 1 



~TT7" + n*°0tr — - - S. 



at 



dG 



dt 



2« w > w >: 



(V, W)) dV 



Po 



dE 



:£, w gocp+((A,W)))dV= I p (((B,W)))dV + f «T,W))dA 

/ JU JdU 



dgo ip 

From this and arbitrariness of W we conclude that [Yavari, et al., 2006] 



dpo . 1 



fdG 

at 



DivP + / 9 B = i o A, 
dE 

2po^ = t, 

og o <p 



(2.24) 

(2.25) 
(2.26) 
(2.27) 
(2.28) 



where P is the first Piola-Kirchhoff stress and r = Jer is the Kirchhoff stress. It is seen that instead of 
conservation of mass we have a balance of mass and the remaining balance laws are unchanged. Note, however, 
that divergence and acceleration both explicitly depend on G, i.e. the time dependency of material metric 
affects the governing balance equations. We will see examples in §2.9. 

2.3 Local Form of Energy Balance. 

Let us now localize the energy balance. First note that 



d „ dE dN dE dG dE T „ dE T 

— E = L V E = 1 : 1 : L V F H : L V K ° tp. 

dt V dN dt dG dt dF v dg VS v 



Note that LvF = because for an arbitrary Z G TyJ5 



L V F=^*(F.Z) = |^Z) = -§Z = 0. 



(2.29) 



(2.30) 



2.4 The Second Law of Thermodynamics and Restrictions on Constitutive Equations 



8 



Using this and also noting that because the background metric is time independent, we have 

d _ dEdN dE dG dE 

d- t E -mlF + dG : ^r + ^ :d ' (2 ' 31) 



where d = ^£vg is the rate of deformation tensor. 4 We know that H = — {^Q,Nyy and [Yavari, et al., 
2006] 

/ ((T,V))dA = / (((DivP,V))+-r:n + T:d)dV, (2.32) 

where Q a b = ^{V a \t, — V(,u), and r is Kirchhoff stress. Thus, from (2.18) and using balances of linear and angular 
momenta we obtain the local form of energy balance as 

dE „ dE dG , 

Po —+BwQ = p a —: — + T:d + Pa R. (2.33) 

In term of the first Piola-Kirchhoff stress, this can be written as 

dE dE dG 

Po— +DivQ = p — :— + P : V V + p i?, (2.34) 

where P : V V = P aA V a \ A . 

Material Metric Evolution. Evolution of material metric is assumed to be given through a kinetic equation 
of the form 5 

dG 

— = T(X,G,F,g)=*(X,G,C). (2.35) 

See Ambrosi and Mollica [2004]; Loret and Simoes [2005]; Fusi, ct al. [2006]; Ambrosi and Guana [2007] for 
some examples written in terms of evolution of F g . We will come back to this problem after first discussing the 
Second Law of Thermodynamics for a growing body. 



2.4 The Second Law of Thermodynamics and Restrictions on Constitutive Equa- 
tions 

In the absence of growth, entropy production inequality in material coordinates has the following form [Coleman and Noll, 
19631 



d 
dt 



p NdV > 



PoR 

e 



dV 



on 



H 

e 



dA, 



(2.36) 



where N = N(X,i) is the material entropy density and = 0(X, t) is the absolute temperature. This is called 
the Clausius-Duhem inequality. 

When material metric is time dependent using balance of mass the Clausius-Duhem inequality is modified 
to read 6 



I ^dV 



dt 



e 



I dA 



N 



dp 
dt 



1 



nPotr — 



2 



<9G 



dt 



dV 



This inequality can be localized to read 



dN Po R 



Pa 



dE dG 
dG : ~dt' 



4 This is the symmetric part of Vv, i.e. the symmetric part of the so-called "velocity gradient". 

5 A simple example that has been extensively studied is the Ricci flow [Hamilton, 1982], which reads 

dG 
~~dt 



2R, 



where R. is the Ricci curvature of the metric G. Note that this flow smooth outs a distorted 3-spherc [Topping, 2006]. 
6 See also Epstein and Maugin [2000]. 
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Note that E = * + N6 and hence 7 

d 1 fdE dtf\ 6 , , 

— N = — )-—(E-m. 2.39) 

dt e \ dt dt J e 2V ' y ' 

Substituting this into (2.38) yields 

Po dE Po d* m ^PoR n . /QV 

e a" " e ■* " Po e* (jE " * } ^ "e" " Dlv { e J + Po ac : aT (2 ' 40) 

Now substituting the local energy balance (2.34) into the above inequality we obtain 

d\P ■ 1 

PiVoV-po^-poNG^-de-Q. (2.41) 

We know that W = *(X, 9, G,F,go tp) and thus 8 

d* d* 3* dG d* , 

Therefore, (2.41) is simplified to read 

/a* \ ■ ( a* \ 1 5* • 

Po ( v ae + N J e+ ( v Po aF _p J :VoV+ e de ' Q + /9o aG :G -°- (2 ' 43) 

Following Coleman and Noll [1963] and Marsden and Hughes [1983] we conclude that 

a* a* 

— = -N and Po^ = P. (2.44) 
and entropy production inequality reduces to po§^ '■ G + ^dO ■ Q < 0. 

Remark. There have been objections in the literature on using the Clausius-Duhcm inequality in continuum 
mechanics [Green and Naghdi, 1977; Marsden and Hughes, 1983]. Next, we show that in growth mechanics 
energy balance and a more general notion of covariance are enough to obtain the restrictions (2.44) on consti- 
tutive equations. In passing we should mention that Green and Naghdi [1991] were able to obtain their entropy 
balance using energy balance and invariance arguments in the case of Euclidean ambient space. What we will 
show next is consistent with their results. 



2.5 Restrictions on Constitutive Equations Using a Thermomechanical Covari- 
ance of Energy Balance 

In this subsection, we follow Marsden and Hughes [1983] and obtain the restrictions (2.44) on the constitutive 
equations using covariance of energy balance with no reference to the entropy production inequality. In the 
case of classical elasticity, Marsden and Hughes [1983] started with the local form of energy balance and postu- 
lated its covariance under simultaneous action of time-dependent spatial diffeomorphisms and time-dependent 
monotonically increasing temperature rescalings. Here we start with the integral form of the energy balance. 



7 Note that * = *(X, 0, G, F, g o tp) and E = £(X, N, G, F, g o tp). 

8 Nota that because is a scalar its time derivative is equal to its Lie derivative along the velocity vector field, and thus when 
the material metric is time independent we have 

d 8W- 9* 9G 9# T „ a* T 9* ■ 9* dG 9* , 
— * = Lv* = © H : 1 : L V F H : L v g = 6 H : 1 : L v g. 

dt ae og at 9F a g 8 9© 9g at a g 8 

We also know that the same time derivative is equal to the covariant derivative of N with respect to velocity vector, i.e. 

d d T „ T 3* A S* 3G 3t „ „ » „ 9* ■ 9* dG 9* „ „ 8*. 9f 3G 9f „„ 

_<I< = — *+V v 4 l = ©H : 1 : V V F-| : V v g = ©H : 1 : V a F = 0H : 1 : VnV, 

dt dt v 96 9G dt 9F v 9g s 9© 9G dt dF m 90 9G dt dF u ' 

where we used the fact that -StW = and V a F = VrjV. See Nishikawa [2002] for a proof. 



2.5 Restrictions on Constitutive Equations Using a Thermomechanical Covariance of Energy Balance 10 

Let us consider spatial diffcomorphisms £ t : S —> S and monotonically increasing temperature rescalings 



s t : m. ■ — > jr. ■ . We assume that at t = to, £ = 1 and ^Ct = We also assume that under these transformations, 
energy balance is invariant, i.e. 



d 

d7 



jf | p' «<B', V')> + R') + p' ^ : ^ + 5; n (V + ± ((V, V'))) | 

+ / (((t',v')} + #Va 

Jaw 



dV 



(2.45) 

Note that £ = ^ + ON. In the new frame (ft = £ o ip and O' = (Q. We assume that E transforms tonsorially, 

(2.46) 



l.C 



E'(X, N', G, F', g') = E(X, N, G, £*F, f.g). 
The same transformation is assumed for free energy density and hence 

d d 9* • 9* • 

s ** . !*(X,Ce,G,£.F, { , g ) = ^6' + JL : G + ±_ : V „V<, 

where (V V) a A = l/ a | A . Therefore, at t = t 

d . a* 

lTr' 



^ = ^ ( e + ,e ) + ||:G + §:(VoV + VoW). 



Thus, at t = to, we can write 
d d 



a<J 



df 



dt dF 



a* 



= —E + — : V W + z — + N 6 + —6' - —6 
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dN' dN 



df 



df 



(2.47) 



(2.48) 



(2.49) 



Energy balance in the new frame at t = to is simplified to read 



dpo 1 
dt 2 



po tr 



A) 



a_G 

a* 



B + ^ ((V, V)) + ((V, W» + \ «W, w» 



^E + ^:VoW + z{ 9 ^ + u]e + ( d ^e'-^e 



dt OF 



ae 



Po (((B'| t=to ,V + W))+i?'| t=to )+p 



dt 

dE dG 
dG 1 ~dt 



dt 



(V + W,A'| t=t 



dV 



S m [E + -{(V,V) 



(V,W)) + -((W,W) 



dy+ / (((T,V + W))+F'| t=t0 )dA (2.50) 

JdU 



Assuming that B' A' = £*(B — A) and subtracting balance of energy (2.18) from (2.45), we obtain 



dpo 1 /5G 
— + oPotr I — ) -S„ 



Si 2' 

a* 



(V,W)) + -((W,W)) 



Po 



u 



dF 



V W + z 



a* 
ae 



n e 



dt dt 



«W,A)) 



dF 



= f po[{{B,W)) + (R'-R)\ t = to ]dV 
Ju 

((Div P, W» + t : (VW) b - (Div Q' - Div Q) t=t(j 



dF, 



(2.51) 



where r : (VW) = t W a i b . Assuming that i? and Q arc transformed such that [Marsdcn and Hughes, 1983] 



dN' 



dt 



dN 



— 9' - R' = C -i-e - i? and Q' = C£*Q 



df 



(2.52) 



2.6 Covariance of the Entropy Production Inequality 



11 



we obtain 



at 



2 PotY 



f dG 
\~dt 



(V, W)) + I ((W, W))) dV - J ((Div P + p B - p A, W)) dV 



: V W-r: (VW)' 



zpo 



9* 



N 9 



dV = 0. 



Note that 



(V W) Q A = g ab [(VW) 
Therefore, arbitrariness of U, W, and z implies that 9 

1 



be 



dpo 

at 



:Po tr 



dG 
~dt 



DivP + / 9 B = p A, 



T T = T 



Po 



dF 



de 



(2.53) 
(2.54) 

(2.55) 

(2.56) 
(2.57) 

(2.58) 
(2.59) 



Thus, we have proven the following proposition. 



Proposition Covariance of energy balance under spatial diffeomorphisms and temperature rescalings gives 
all the balance laws and the constitutive restrictions imposed by the Clausius-Duhem inequality. 



2.6 Covariance of the Entropy Production Inequality 

In this subsection we study the consequences of covariance of the Clausius-Duhem inequality. Again, let us 
consider the diffeomorphisms £t : S — > S and monotonically increasing temperature rescalings Q '■ K + - * 
Let us postulate that the entropy production inequality is invariant under the simultaneous action of these two 
transformations, i.e. 



d 
rff 



u 



PoN'dV > 



u 



p^R> 
6' 



dV 



m 



H' 
W 



dA 



N' 



u 



dt 



1 



Potrl 



fdG 



dV 



u 



,dE> dG 



Note that (2.52) implies that 



(2.60) 
(2.61) 



dN' R' dU H' H 

-df-^ = ^- R and e^=e- 

It can easily be shown that the inequality (2.60) is identical to (2.37), i.e. assuming the transformations (2.52), 
entropy production inequality is trivially covariant. 



'In a previous footnote it was shown that 



Thus 



d 9# A 9* dG 9* „ , r 9* • d<b dG 9* „ 

— * = H : 1 : V V = 6 H : 1 : £ v g- 

dt 90 dG dt dF dB dG dt 9g & 

9* 

™ : V V = — : £ v g. 
9F 9g 



This holds for an arbitray change of frame ^ : 5 -> 5 as well, i.e. at t = to: 



9* 
9F 



V (V + W) 



9<B 
9g 



£(V+W)g- 



Hence, for arbitrary W 



Noting that [Marsden and Hughes, 1983] 
Doyle-Ericksen formula. 



dE_ 

0(p 



d$> „ „ T dV „ 
— : V W = — : £ w g. 
9F 9g 

, this means that po §p = P is equivalent to 2po 



i.e. the 
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2.7 Principle of Maximum Entropy Production 



In this subsection we use the so-called maximum entropy production principle to obtain a kinetic equation for 
G. This principle states that a non-equilibrium system with some possible constraints evolves in such a way to 
maximize its entropy production [Ziegler, 1983; Rajagopal, K. R. and Srinivasa, A. R., 2004a]. This principle 
has found applications in many different fields of science. For a recent review see Martyushev and Seleznev 
[2006]. This principle has recently been used in growth mecahncis for obtaining kinetic equations for "growth 
velocity gradient" [Loret and Simoes, 2005; Ambrosi and Guana, 2007; Fusi, et al., 2006]. Here, we use it in our 
geometric framework. 

For a growing body, entropy production in a subbody 1A C B is defined as 



R 



dt 



'e 



Q 



Po-jt - Pott + Div -pr - : -57- 



e 



- / n 

dU 

dE dG 



dp 



dt 



dG 



dt 



dV - 



dE dG 



'dG dt 



A 



dV= I -dV, 



(2.62) 



where A 
we have 



(-) 



Pof -Pof+Div 



A 



n dE . dG 
POOG ■ dt 



dG ' dt 

is the rate of entropy production. Using energy balance, 



-Po 



d$> 

w 



1 



P: V V--de-Q-9N. 



Note that 

Thus, we can write 



dtf . a* dG „ „ Tr 

-p ^^eN-p -:--P:V V. 

<9* • 1 



(2.63) 
(2.64) 

(2.65) 
(2.66) 

X+l-dG- ^ 

Note that part of entropy production rate is constitutively given, i.e. A = A(0, G, G, F, g) — gd8 ■ Q. As the 
simplest example let us assume that 



We now maximize A with respect to G under the constraint (2.65). Let us define 

/ a* • 1 

$ = A + A I A + Pa — : G + -dG • Q 

where A is a Lagrange multiplier. Maximizing $ with respect to G gives 



dA 
9G~ 



A 



-Po- 



A = /? tr G 2 - id6 • Q = (3G ab G C dG ac G bd - ^d9 • Q. 



Thus 



Or 



dA 

ddZ 

2/3 dG 



2{3G AC G BD G CD . 



A 



Using (2.70) and (2.68) we can write 

A = 



A 



9* 

~2/?(A + l) P °dG' 



4/3(A + l) 2 ™<9G ■ dG 
At the same time using (2.70) and (2.65) we have 



^de • q. 



A = 



A 



a* a* 1 



2/3(A + lj Po 'dG ' dG 
Looking at (2.71) and (2.72) we see that A = —2 and hence 

r* 1 9 * 
= -j3 P °dG- 



e deQ - 



(2.68) 
(2.69) 
(2.70) 

(2.71) 

(2.72) 

(2.73) 
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2.8 Isotropic Growth 

For isotropic growth, material metric has the following time dependent form: 

G(X, t) = e 2O(x ^ ) G (X), (2.74) 
i.e. a family of conformal material metrics model the growth. Thus 

—^— = 2-G(X,t). (2.75) 
Therefore, balance of mass is simplified to read 



^ + p (X, i )«e^)trG (X)^ m (X, t ). 
Given G = G(X,t), one has the following relation between volume elements at to and t: 



(2.76) 



Or 

dV(X, t) = e Nn ^dV (X), (2.78) 
where N = dimS. Mass form has the following representation 

m(X,t) =po(X,*)dV(X,t) = e Jvn(x '* ) m (X). (2.79) 



Note that rh(X, t) = iV^Pm(X, i). Mass of a subbody U C B will have the following time-dependent form 



U t (U) = [ m(X,t). 
Ju 



(2.80) 



Hence 

^M t (W) = jN^m(X,t). (2.81) 

In isotropic growth there is no change in shape due to addition or removal of mass. In the decomposition of 
deformation gradient one has 

F = F e F 9 and F 9 =g(X,f)I, (2.82) 

where I is the identity map and g is a scalar field. We will discuss this decomposition in more detail is §3. In 
the sequel we will obtain those isotropic growth distributions that are stress free. But let us first look at some 
simple examples of material metric evolution. 

2.9 Examples of Bulk Growth 

In this subsection we look at three examples of bulk growth and show how analytical solutions for residual 
stresses can be generated for both isotropic and non-isotropic growth. 

Example 1 (Isotropic Growth of a Neo-Hookean Annulus): Let us consider a two-dimensional, in- 
compressible neo-Hookean material in a flat two-dimensional spatial manifold. The free energy density of a 
nco-Hookcan material in two dimensions has the form 

# = ¥(X,C)=M(trC-2), (2.83) 

where C is the Cauchy-Green tensor, or equivalcntly, the pull-back of the spatial metric, Cab = F a AF b B9ab, 
and fi is a material constant. We assume that this form holds for a growing isotropic material. In components 



* = n (F a A F b B9ab G AB - 2) 



(2.84) 
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14 



The "2" is of no particular significance: when the material metric is fixed, it simply shifts the free energy by 
a constant. When the material metric changes its contribution to the free energy is proportional to the time- 
dependent material volume, which, for a given growth distribution, is independent of the spatial configuration. 
We ignore this term, and use 4" = ^(X, C) = fi tr C as our definition of the free energy. 

Let us assume that initially the material has a flat annular shape R\ < R < R2 without any stresses. We 
would like to calculate the stresses that occur in the new equilibrium configuration after a rotationally symmetric 
growth, CI = Q(R,t). In polar coordinates, the spatial metric and its inverse read 

and thus detg = r 2 . The only nonzero connection coefficients are: 7^ = — r, 7^ = 7^, = 1/r. For the 
rotationally symmetric time-dependent material metric we have 

r _ ( G RR G Re \ _ 2U(fi,t) ( 1 \ x _ / G RR G R& \ _ 2 n(fi,t) ( 1 \ , , 

G -{ Ge R Gee )~ e [OR 2 )' G ~{ G eR G ee ) e [o l/R 2 ) ' (2 ' 86) 

and thus, det G = R 2 e iQ( - R ^. The following nonzero connection coefficients are needed in the balance of linear 
momentum: 

T R R = Q'(R, t), rg e =-R- R 2 Q'(R, t), ri e - r| fl = 1/R + n'(R, t), (2.87) 
where H.'(R,t) = Given Q = Q(R,t), we are looking for solutions of the form 

<p(R,Q) = (r,6) = (r(R,t),Q). (2.88) 



,10 



Thus 



This gives the Jacobian as 



p=( ^ M, p-i = ( Vr'CR,*) o j (2 89) 



Incompressibility dictates that 



rr' = Re™^. (2.91) 
This differential equation has the following solution 



R 

r 2 (R,t) = r((R,t) + I 2^e 2n ^' t) dt (2.92) 

JRi 

Note that r\ (R) is not known a priori and will be obtained after imposing the traction boundary conditions at 
r% and r 2 - In incompressible elasticity, P aA is replaced by P aA — Jp(F~ 1 )~ A i,g ab , where p is an unknown scalar 
field (pressure) that will be determined using the constraint J = 1 [Marsdcn and Hughes, 1983], i.e. 

P aA = 2 t iF a B G AB -p(R)(F- 1 ) A b g ab . (2.93) 

Therefore, using (2.91), we obtain the nonzero stress components as 

p- R = 2 J^--p(R)L e -^W and = ^e-anCH,*) _ pW (2 94) 

r R Rr 

where p(R) is an unknown pressure. 



10 If one does not consider the intrinsic metric and instead uses the standard metric of the Euclidean space, F has the following 
representation 

F = 



r'(R,t) \ / e "C«.t) o 



Thus, detF e = and hence J e = 1 is equivalent to (2.91) as expected. 
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Balance of linear momentum in components reads 

r A 1 AB + P lbc F A = 0. 



paA\ u i j~>A t nbA_,a 



(2.95) 



For the radial direction, a ~ r, we have 



rA 



flP 

prAl KJ± i tA jjrB , TjbA„,r j^c 



dX A 

QprR 

dR 

QprR 

dR 



n B P^+P DA Y bc F C A 



RR 1 0R ) 
1 

R 



2n'(R,t) P rH -rP 



?0O 



This gives 



p'(R) 



2fiR 



,2n(R,t) 



2(1 + RQ') 



Assuming that p{Ri) = 0, we obtain 



p(R) 



r"(0' 



2(1 + 



E>2 „2 
£_ 2f!(Ji,i) _ ' -2Q(fl,t) 

r 2 i? 2 



2fi(£,t) _ ^(0 --2n(C,t) 



For a — 8, balance of momentum (2.105) gives 



P* A U = 



96 



r A eP ,e + jrf^^ + P oe lSgeF s e = (rfe + r g e) p *e = q 



(2.96) 



(2.97) 



(2.98) 



(2.99) 



i.e. this equilibrium equation is trivially satisfied. Note that tr (^f) = 40' and hence balance of mass reads 

d Po (R,t) 



at 



2Cl'(R,t)p (R,t) = S m (R,t). 



(2.100) 



This differential equation can be easily solved for mass density. 

Note that if one considerers a cylinder with kinematics assumptions r — r{R, t),6 — 6, z — kZ for a constant 
k, the residual stresses will be very similar to what was just calculated. 

Example 2 (Anisotropic Growth of a Neo-Hookean Annulus): Let us consider an anisotropic growth 
represented by the following material metric: 



G 



3 20(fl,t) 







£>2g-2fi(H,t) I i 



G 1 = 



e -2n(R.t) q 

1/R 2 e 2n ^ J ' 



(2.101) 



and thus, det G = R 2 . The following nonzero connection coefficients are needed in the balance of linear 
momentum: 



r* R = tf(R,t), T R e =Re- An ^[n'(R,t)-l], T% @ =T% R = l/R-il'(R,t). 



(2.102) 



Given £1 = fl(R,t), we are looking for solutions of the form tp(R,Q) = (r,6) = (r(R,t),Q). Thus, F and F 1 
have the forms given in (2.89)and this gives the Jacobian as J = L R -. Incompressibility dictates that rr' = R. 11 
This simple differential equation has the following solution 



r{R, t) = y/R 2 + C(t) = ^R 2 -Rl + rl 



(2.103) 



'In the classical formulation 



r'(R,t) 

r(R,t) 
R 







and F„ 









-f!(fl,t) 



Thus, J e = det F e = 1 would lead to the same incompressibility constraint rr' = R. 
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Note that T\ (R, t) is not known a priori and will be obtained after imposing the traction boundary conditions 
at ri and r^. Now, we get the nonzero stress components as 



P rR = 2f,e- 2n ^ 



R 



r(R,t) 



■p{R,t) 



R 



and 



P ee = 2^ 



2Q(R,t) 



p(R,t) 
r 2 {R,ty 



(2.104) 



where p{R, t) is an unknown pressure. 

Balance of linear momentum in components reads 



paAl 



For the radial direction, a = r, we have 



jrA\ 



This gives 



P '(R,t) 



2fiR 



Assuming that p(Ri, t) = 0, we obtain 
p(R,t) = 



2tf 



Note that r 2 = R 2 + C and thus 



, e + c' 



QpaA 

dX A 



QprA 

dX A 

gprR 

dR 

gprR 

dR 



-2Q(R,t) 



-2fi(f,t) 



-2Q(£,t) 



T%P aa + P bA ^ c F c A = 0. 



(2.105) 



Ti B P rB + P b HcF c A 



RR 



0R. 



P 



rR 



P 



± p rR _ rP oe = 
R 



2 - 2RQ'(R,t) 



2-2£fi'(£,i)- 



R 2 ' 



4Q(R,t) 



R 2 



r (6 „4n(g,t) 



r 2 (£) 



e + c 



(2.106) 
(2.107) 

(2.108) 
(2.109) 



Assuming that p(i?2, t) = 0, C(t) can be calculated using the above equation. 
For a = 9, balance of momentum (2.105) gives 



P eA U 



dp ee 

80 



Ti e P s ° + P 9R -i e rr F r R + P^ ee F° e = (r« e + rg e ) P 



0. 



(2.110) 



i.e. this equilibrium equation is trivially satisfied. It is seen that a growth that results in only change in shape 
and no change in volume can still result in residual stresses. Note that tr (^) 
reads 

dp (R,t) 



at 



S m (R, t). 



and hence balance of mass 
(2.111) 



Example 3 (Spherical Growth of a Neo-Hookean Hollow Sphere): Let us consider a hollow sphere 
with inner and outer radii Ri and R initially in a coordinate system (i?, 0,$). Let us denote the spatial 
coordinates by (r, 9, </>). The spatial metric has the following form 



1 

Or 2 
r 2 sin 2 <; 

with the nonzero connection coefficients jg g = —T^ r f±± = — r sin 2 (j>, 7^ = lg = ^l r ilt 
isotropic growth of the hollow sphere we consider the following material metric 

10 
= e 2 ^*'*) I R 2 

R 2 sin 2 $ 



(2.112) 
1/r. For 

(2.113) 
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The nonzero connection coefficients are 



r* fl = a', rg e = -r- r*. = -(r + R 2 n<) sin 2 $, r« e = rg^ = r£ $ = rf fl = - + n'. (2.114) 

Under this symmetric change of material metric (growth) we look for solutions of the form r = r(R,t), 8 = 
6, = Thus 12 

r'(R) 

1 ) , (2.115) 
1 



and hence 



Incompressibility gives us 



The only nonzero stresses are 



R 2 



r 3 (R)=r 3 (R) + [ 3£ 2 e 3n («'<^ 



yrR 



« ^ n r 



-3f2 pSB 



2fi 



2n _ _P_ p<>* 

O 7 



_ e -2Q 



i? 2 sin 2 $ 



r 2 sin 2 $ 



(2.116) 
(2.117) 

(2.118) 



Again, the equilibrium equations for P ee and P^ are trivially satisfied. The only nontrivial equilibrium 
equation reads 



P tA \a 



dP 



rA 



dX A 

QprR 

dR 



r% B P™ + P , ~fl c F c A 

frgfl + r| fl + r| fl ) p rR + p ee Y ee F s e + p** 7 j^% 



^- + 30' ) P rR - rP ee - r sin 2 <S>P^ = 0. 



This gives 



p'{R,t) 



dR \R 



4flR AQ(Jj^) 



I + Q'(R, t) - ^e 3n ^ - -^e~ 3n ^ 



Assuming that p(R±) = 0, we obtain 

P (R,t) = 



R „4f2(g,t) 



l r 4 (0 



r 3 (£) 



Note that tr (^) = 6f2' and hence balance of mass reads 



dpp(R,t) 
dt 



3n , {R,t) Po (R,t) = S m (R,t). 



(2.119) 



(2.120) 



(2.121) 



(2.122) 



This differential equation can be easily solved for mass density. 



12 In the classical formulation 



/ r'(R,t) \ 

R 



and F„ 



o 

V o o tim J 

Hence, J e = det F e = 1 would lead to the same incompressibility constraint. See Chen and Hogcr [2000] for more details. 



e n( - R >V 
e n («>*) 
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2.10 Visualizing Material Manifolds with Evolving Metrics 

In our geometric theory, we model growth in a fixed material manifold B. We can visualize the evolution of 
G(t) by embedding B in some material ambient space X with a fixed metric H. For us this larger space would 
be the Euclidean space with its standard metric. Consider a one-parameter family of isometric embeddings 
it : B X, i.e. t£H = G(t). For the sake of simplicity, let us restrict ourselves to rotationally symmetric 
metrics, i.e. we look at metrics of the form 

G _(M\R,t) \ 

N 2 (R,t) ) ' [2 - U6) 

in some coordinate patch (R, 0), i.e. metric has the form: M 2 (R,t)dR 2 + N 2 (R,t)d® 2 , where t is time. Note 
that M and N are independent of O. We now look for solutions in the set of surfaces of revolution. Let us 
consider a time-dependent curve 7(s,t) = (p(s,t), £(s, t)) in the plane. The surface obtained from this curve by 
revolution about z-axis has the following parametric representation: 

$(s,9,i) = (p(s,i) cos e,/o(s,t) sin e,£(s,i)). (2.124) 

The induced Euclidean metric is [Peterson, 1997]: 

(dX 2 + dY 2 + dZ 2 ) = (p 2 (s, t) + i 2 (s, tf) ds 2 + p(s, t) 2 dB 2 , (2.125) 

where a superimposed dot means differentiation with respect to s. Given M(R,t) 2 dR 2 + N(R, t) 2 d<d 2 , let us 
assume that p(s,t) = N(s,t) and hence 



£(s,t) = ^M 2 ( Sl t)-N 2 (s,t). (2.126) 

Therefore 



£(*,*)=/ yM 2 {£,t)-N 2 (e,t)d£. (2.127) 

J so 

Of course, a solution may not exist. This happens when M 2 < N 2 . This is not surprising as not every 
rotationally symmetric metric arises from a surface of revolution. In the following we consider an initially 
stress-free annulus under different rotationally symmetric growth distributions. 

Example 1: Consider isotropic growth, i.e. 

G= " (R) , M = e Q , N = Re n . (2.128) 



R 2 e 2n ^ , 

Hence M 2 — N 2 = -Rn'e 2n (Rfl' + 2). Let us look at two cases: 

i) Q{R) = —R: We have M 2 —TV 2 = e- 2R (2R-R 2 ), which for < R < 2 gives the material manifold shown 
in Fig. 2.2(left). 

ii) n(R) = -R 2 : We have M 2 - N 2 = 4R 2 e - 2R2 (l - R 2 ), which for < R < 1 gives the material manifold 
shown in Fig. 2.2(right). 



Example 2: We look at anisotropic metric evolutions represented by 

R 2 e 2Yl ^ 

We look at two cases: 



G= n , M = e u , N = Re n . (2.129) 



i) n(R) = cos 2 R and U(R) = 0: We have M 2 - N 2 = e 2 cos2 R - 1 > 0. The material manifold shown in Fig. 
2.3(left). 

ii) Q(R) = and U(R) = - \nR 2 : We have M 2 - N 2 = 1 - which for R > 1 gives the material manifold 
shown in Fig. 2.3(right). 
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Figure 2.2: Visualization of the material manifolds of two isotropic growth distribution of an annulus as embeddings in R 3 . Left: 
fl(R) = -R. Right: Q(R) = -R 2 . 





Figure 2.3: Visualization of the material manifolds of two anisotropic growth distribution of an annulus as embeddings inR 3 . Left: 
= cos 2 R, U(R) = 0. Right: Q(R) = 0, n(_R) = - \nR 2 . 



2.11 Stress-Free Isotropic Growth 

In the context of growth mechanics, Takamizawa and Matsuda [1990] realized that having a stress-free config- 
uration is equivalent to vanishing of Riemann's curvature tensor, although they did not present any detailed 
calculations. In this subsection we study this problem in detail and obtain stress-free isotropic growth distribu- 
tions in both two and three dimensions. 

Let us first review some basic concepts in Riemannian geometry. For 7r : E — > S a vector bundle over a 
manifold S, £(S) the space of smooth sections of E, and X(S) the space of vector fields on S, a connection on 
E is a map V : X(S) x £{S) -> £(S) such that V /, f u f 2 G C°°(S), V a x , a 2 £ M 

a) V/ lXl +/ 2 x 2 Y = / x V Xl Y + / 2 Vx 2 Y, (2.130) 

b) Vx(aiY 1 +a 2 Y 2 ) = a 1 Vx(Y 1 )+a 2 Vx(Y 2 ), (2.131) 

c) Vx(/Y)=/V X Y+(X/)Y. (2.132) 

A linear connection on S is a connection on TS, i.e., V : X(S) x X(S) — * X(S). In a local chart {x 1 } 

^dA=lpk, (2.133) 

where 7^ are Christoffel symbols of the connection and di = . A linear connection is said to be compatible 
with the metric of the manifold if 

V x «Y, Z» = ((V X Y, Z)) + ((Y, V X Z)> . (2.134) 
One can show that V is compatible with g if and only if Vg = 0. Torsion of a connection is defined as 

T(X,Y) = VxY-V Y X- [X,Y], (2.135) 

where 



[X,Y](F) = X(Y( J F)) - Y(X(F)) V F e C°°(S), 



(2.136) 
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is the commutator of X and Y. V is symmetric if it is torsion- free, i.e. VxY — VyX = [X, Y]. According to 
the Fundamental Lemma of Riemannian Geometry [Lee, 1997] on any Riemannian manifold (S, g) there is a 
unique linear connection V, the Levi-Civita connection, that is compatible with g and is torsion-free with the 
following Christoffcl symbols 

^~2 9 {^ + d^-^)- (2 ' 137) 

Curvature tensor TZ of a Riemannian manifold (S, g) is a f^J -tensor H : T£S x T X S x T X S x T X S — > R 
defined as 

92.(a,wi,W2,w 3 ) - a (V Wl V W2 w 3 - V W2 V wl w 3 - V [wi , W2 ]W 3 ) (2.138) 
for a € T*S, w 1 ,w 2 ,w 3 G T x 5. In a coordinate chart {x a } 

_ dlL d lbc , ni a _,e „,a _,e i o i gn\ 

l< bed = ~ -q^J + Icelbd - Idelbc- (2.139) 

Note that for an arbitrary vector field w 

™V - w a \ cb = n a bc dw d + T d cb w a ld . (2.140) 

An n-dimcnsional Riemannian manifold is flat if it is isometric to Euclidean space. A Riemannian manifold is 
flat if and only if its curvature tensor vanishes [Lee, 1997; Spivak, 1999; Bcrgcr, 2003]. Ricci curvature is defined 
as 

Rab = n C a cb. (2.141) 

The trace of Ricci curvature is called scalar curvature: 

R = R ab g ab . (2.142) 

In dimensions two and three Ricci curvature algebraically determines the entire curvature tensor. In dimen- 
sion three [Hamilton, 1982]: 

Ttabcd = gacRbd ~ g a dRbc ~ gbcRad + gbdRac ~ ^ ^ (ffocSfed — 5ad.9bc) ■ (2.143) 

In dimension two R ab = R.9ah, and hence scalar curvature completely characterizes the curvature tensor and 
is twice the Gauss curvature. Note that any one-dimensional metric is flat. In the following we obtain the 
stress-free growth distributions in dimensions two and three. 



i) The two-dimensional case. Consider a two-dimensional shell restricted to live on a flat planar surface 
between two rigid planes. We assume that with no external or body forces, initially the shell is stress-free. Can 
one find the growth distributions that will result in equilibrium configurations with zero stress? Uniform growth 
will obviously result in uniform expansion/contracction, and hence no stress. Are there other isotropic growth 
distributions with this property? 

The spatial distances between material points are measured by the ambient space metric (the "spatial 
metric"), which is Euclidean. A given growth distribution will result in a change in the material metric. A 
configuration will be stress- free if there is no "stretch" in the material, i.e., if the material distance between two 
points is the same as the spatial distance. This can happen only if the two metric tensors (spatial and material) 
give the same distance measurements between nearby material points, i.e. if they are isometric. As the spatial 
metric is assumed to be Euclidean, this means that the material metric, after the change due to a given growth 
distribution, must be Euclidean. 

Riemann defined the curvature tensor of the metric and proved that a metric is flat, i.e., it can be brought 
into the Euclidean form Su locally by a coordinate transformation, if and only if its curvature tensor is zero 
[Lee, 1997; Spivak, 1999; Berger, 2003]. It turns out that in dimension two, a weaker requirement is sufficient 
[Bcrgcr, 2003]: a metric is flat if and only if its scalar curvature (the Ricci scalar) is zero. Let us now apply this 
condition to a two-dimensional metric that is obtained from a non-uniform growth distribution on an initially 
stress-free, planar shell, i.e., Gjj = e Sij. The Ricci scalar for a metric of this form is given by [Wald, 1984] 



R = -2e~ 2a V 2 ft. 



(2.144) 
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Thus, R = requires V 2 i7 = 0, i.e., fl has to be a harmonic function. Note that here V 2 is the spatial Laplacian. 
Growth is a slow process compared to elastic deformations and therefore time can be treated as a parameter 
and hence incrtial effects can be ignored. Hence, time in Q is treated as a parameter. 

It is worth emphasizing the distinction between local and global flatness, and the implications for stress-free 
growth distributions. Although the surface of a right circular cylinder in three dimensions looks curved, it 
is locally, intrinsically flat. For any given point on the cylinder, one can find a finite-sized region containing 
the point, and a single- valued coordinate patch on this region, for which the metric has the Euclidean form. 
Physically, this means that for any given point, we can cut some finite-sized piece containing the point, and can 
lay the piece on a flat plane, without stretching it. The surface of a sphere in three dimensions, on the other 
hand, is intrinsically curved; it is impossible to make any finite-sized piece of the sphere, no matter how small, 
to lie on a flat plane without stretching it. The curvature condition R = (or V 2 fi = 0) is local. Making a full 
cylinder to lie in a plane nicely (i.e., without tearing, folding, or stretching it) being impossible is due to the 
global topology of the cylinder; local restrictions on curvature cannot constrain the global properties sufficiently. 

Let us specialize to the case where Q depends only on the radial coordinate R of an initially flat annular 
piece of a material, Rq < R < Ri. The flatness condition gives 



V 2 n 



RdR 



(2.145) 



Solving this gives 



3 2n = £(t)R 2,lit) 



(2.146) 



where £ > and rj are time-dependent constants. The metric rescaling (2.146), with the proper identifications, 
is describing an annular piece from a conical surface, with deficit angle £ = 27r(l — l/|c|), where c = 
[Ozakin and Yavari, 2009] . Now, one can show that it is not possible to make such a conical surface lie on 
the plane without tearing, stretching, or folding it. Thus, starting with an annular shell between two rigid 
planes, a growth distribution of the form (2.146) will indeed result in stresses, although the related material 
metric is intrinsically flat (see Fig. 2. 4. a). However, if the material consists only of a simply-connected piece 
of the annulus (say, R\ < R < R2, < 9i < O < 62 < 27r), the growth distribution (2.146) will just cause a 
stress-free expansion of the material, between the two rigid planes. See Fig. 2.4.b. 






a) b) 

Figure 2.4: a) Zero-stress growth of an annulus to a cone, b) Zero-stress growth of a simply- connected piece of an annulus. 



A remark on conformally flat manifolds and growth mechanics. A Riemannian manifold (B, G) is 
conformally flat if there exists a smooth map / : B — > M such that G = f8, where 8 is the Euclidean metric. In 
isothermal coordinates the conformally flat Riemannian metric has the following local form 

G = /(X)(rfA 2 + ... + dA 2 ). (2.147) 

It is known that [Bcrgcr, 2003] any two-dimensional Riemannian manifold is conformally flat and the map / 
is unique. A corollary of this theorem in our theory of growth mechanics is that given any smooth curved 2D 
stress-free solid, there exists a unique growth distribution such that in the new (grown) configuration, the 2D 
solid is flat and still stress free. Equivalently, starting from a stress free flat sheet, it is always possible to deform 
it to any smooth curved shape by growth without imposing any residual stresses. 
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ii) The three-dimensional case. Let us next consider the three-dimensional case. In three dimensions, a 
vanishing Ricci scalar is not sufficient to guarantee local flatness. However, a three-dimensional metric is flat if 

and only if its Ricci tensor vanishes [Berger, 2003]. The Ricci tensor Ru of the metric Gu = e 2 ° Gu is given 

o o 

in terms of the Ricci tensor Ru of Gu by the following relation [Wald, 1984] 

Ru =Ru -(n — 2)V/VjO- GuG V^V^O + ( n - 2)V/OVjO - (» - 2) GuG V K nv L 0, (2.148) 

o o 

where n is the dimensionality. Now, once again, assume that the initial metric Gu= Su, Ru= 0, and n = 3, 
and replace the covariant derivatives with partial derivatives. This gives 

Ru = -didjd - 5 u & KL d K d L 9, + djOdjCl - 5u5 KL d K nd L n = 0. (2.149) 

This gives the following system of nonlinear partial differential equations in terms of O: 

0,12 = 0,10,2, (2.150) 

0, 13 = 0, 1 0, 3 , (2.151) 

0,23 = 0, 2 0,3, (2.152) 

0,n + V 2 + 2 2 + 2 3 = 0, (2.153) 

0,22 + V 2 + O 2 ! + 2 3 = 0, (2.154) 

0,33 + V 2 + O 2 ! + 2 2 = 0. (2.155) 

This system of nonlinear equations were solved in [Ozakin and Yavari, 2009] . The general solution is 

n{X\ X 2 , X 3 , t) = -In (co(t) [(X 1 ) 2 + (X 2 ) 2 + (X 3 ) 2 ] + Cl p' + c 2 (t)X 2 + c 3 (t)X 3 + c 4 (t)} . (2.156) 

In a special case if c\ = c-2 = C3 = C4 = 0, we have 

n(X\X 2 ,X 3 ) = -In (c R 2 ) , (2.157) 



where R = y/ (X 1 ) 2 + (X 2 ) 2 + (X 3 ) 2 . In order to understand what this solution represents physically, let us 
write the metric in polar coordinates. 



dS 2 = e 2n [dR 2 + R 2 (dG 2 + sin 2 9d$ 2 )] = [dR 2 + R 2 (dQ 2 + sin 2 9d$ 2 )] . (2.158) 



Now let us define 



R=± (2.159) 



In terms of i?, the metric becomes 

dS 2 = dR 2 + R 2 (dQ 2 + sin 2 6# 2 ), (2.160) 

which is precisely the flat Euclidean metric in three dimensions. Thus, after the growth, the metric is still flat, 
but the radial coordinate in which it is manifestly so is related to the old radial coordinate by (2.159) (up to 
a simple shift of origin). This means that, particles at the two radii Ri < R2 move to the new radii i?i > R2, 
after the growth, i.e., the material gets "inverted". This may not be possible for a solid ball without tearing it 
apart, but it is perfectly possible for a piece from such a ball. 

If only C4 is nonzero, we recover the trivial uniform growth. If only c\ is nonzero and assuming that the 
initial material metric is Euclidean for the half space X 1 > 0, we have 

Gu = (Yrj^"' ( 2 - 161 ) 
where A = l/(ci) 2 . This shows that the material manifold is conformal to the Poincare half space. 
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2.12 Lagrangian Field Theory of Growing Bodies 

In the Lagrangian formulation of nonlinear elasticity, one assumes the existence of a Lagrangian density 

C = C(X,t,G,ip,ip,F,g). (2.162) 
Lagrangian is defined in the reference configuration as 

L= f C(X,t,G(X),ip(X),<p(X),F(X),g(ip(X)))dV(X). (2.163) 

JB 

In the case of a growing continuum material metric will be a dynamical variable too. Thus, for growth of an 
elastic body we assume the existence of a Lagrangian density C = C (X, t, G, ip, ip, F, g) and write the Lagrangian 

as 

L = J £(x,t,G(X,t) ) ^(X,i),< / j(X,t),F(X,t),go^(X,t))dF(X) J (2.164) 



where cfV(X) = v 7 det G dX 1 A ... A dX n = vdetG <iX. Having the Lagrangian, action is defined as 

S = [ Ldt (2.165) 



J to 

and Hamilton's Principle of least action states that 

SS = dS ■ (Sip, SG) = 0. (2.166) 

The problem with this formulation is that it assumes that the solid is a conservative system. This is obviously 
not correct here as growth is a dissipative process, in general. There have been recent works on Lagrangian 
formulation of dissipative systems. One idea is to use fractional derivatives and assume that Lagrangian is a 
function of some non-integer time derivatives of generalized coordinates [Riewe, 1997]. It is not clear how one 
can use this idea for a general field theory and even if successful how useful that theory will be. Another way of 
considering dissipation in Lagrangian mechanics is to use a Rayleigh dissipation function [Marsden and Ratiu, 
2003]. 

Assume that there exists a Rayleigh dissipation function TZ = TZ(G, G). For a continuum with dissipative 
forces F, the Lagrange-d'Alembert Principle states that [Marsden and Ratiu, 2003] 

S [ [ CdVdt + [ [ F • 5(p dt = 0. (2.167) 

Jto JB Jto JB 

Assuming the existence of a dissipation potential TZ for a growing body, the two dissipative forces are represented 
as 

F = -^ and F G = -^2. (2.168) 
dip dG 

In this case, Lagrange-d'Alembert Principle states that 



I' [ C(X,t, G,(p,<p,F,goip)dVdt+ I' { (F • Sip + F G • SG) dVdt = 0. 

Jt JB Jt JB 



(2.169) 



For the sake of simplicity, let us consider the two variations separately. 
Case 1: If only deformation mapping is varied, one has 

SS = dS ■ (Sip, 0) = 0. (2.170) 

This can be simplified to read [Yavari, et al., 2006] 

dC d dC ( dC \ dC nc b n dC b dTZ ,„,„,. 

W a " dtW " \dF^l ) lA " dpr A F Alac + 2 d^ d 9bdlac = w ( } 
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Or 

P a A \A + + (F c AP b A - Jo cd g hd ) t* c = Po9ab A b + |p. (2.172) 

Case 2: In material representation of classical nonlinear elasticity density is not a dynamical variable but rather 
it is a parameter appearing in the Lagrangian. It is through a reduction process (material to spatial) that the 
density ends up satisfying the continuity or advection equation (see Holm, et al. [1998] for more details). Here, 
we should note that unlike classical nonlinear elasticity, mass density varies, in general, when material metric 
changes. In other words, 5po and 5G arc related through the nonholonomic constrant of mass balance. For a 
similar discussion on Lagrangian formulation of fluid mechanics in Eulerian (spatial) coordinates see Brcthert 
[1970]. We know that by definition of S m 

d 
dt 



Po (X,t)dV= / S m (X,t)dV= / S m (X,t)dV, (2.173) 



where Sm (X,t) is mass source in the initial material manifold with volume element dy . Note that S is the 
quantity that can be given physically. Now balance of mass can be rewritten as 



Po (X,t)dV= / p (X,t )dV + / / S m (X,t)dV dr. (2.174) 

JU J t JU 

o 

Now for a fixed Sm, let us consider mass density and material metric variation fields po(X, i; e) and G(X, t; e). 
For an arbitrary e the above integral mass balance reads: 

/ Po (X,t;e)dV e = / Po (X,t ;e)dV e + f f S m (X,r)dV dr. (2.175) 

JU JU J t JU 

Let us take derivatives with respect to e of both sides, evaluate them at e = and note that all variations vanish 
at t = to- This gives us 

$Po + ^Po tr{SG) ) dV = 0. (2.176) 



As U is arbitrary we obtain 

5po + ^p tr(SG) = 0. (2.177) 
We can write C = PqC, where C is Lagrangian density per unit mass. Hence 



ti 

' dVdt. (2.178) 



CdVdt = 6 / poCdVdt = 

to JB Jto JB Jt JB 



Spa + \po tr(<5G) 



Thus, using (2.177) 

5 1 [ CdVdt = I [ SCdVdt, (2.179) 

Jto JB Jto JB 

where in 5C mass density is assumed to be fixed. Now substituting this in the Lagrange-d'Alembert Principle 
and assuming that G varies while if is fixed, we obtain 

U - ^] : SGdVdt = 0. (2.180) 
dG dGJ 

Therefore, the corresponding Eulcr-Lagrange equations are 

f£ = ^. (2.181) 
3G 3G 

Note that this is very similar to what we obtained using the principle of maximum entropy production in §2.6. 
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Example. Assuming that Tl(G, G) is a quadratic function, i.e. TZ(G, G) = wtrG 2 = luG A mG B nG ab G mn , 
we have 

= 2loG M nG am G bn . (2.182) 

OGab 

Thus 



1 dC 

W = 

Note that £ = T - V, where 



G -^dG- ^ 



T=\po((V,V)) and V = Po E + V B , (2.184) 
where V_b is the potential of body forces. Therefore, the evolution equation (2.183) reads 

G» - -i*J§. (2.185) 



Note that E = * + N6 and hence 



Thus 



dE /ov ay oe\ 9* 

9G = Ug + M5G +N 9G-9G- (2 ' 186) 



G« = -^o|§, (2.187) 



which is identical to (2.73) if we choose u) = h(3. 

3 Connection between F = F e F 5 and the Geometric Theory 

In the literature of growth mechanics the idea of multiplicative decomposition of deformation gradient into elastic 
and growth parts is usually attributed to Rodriguez, et al. [1994], although it can be seen in several earlier works 
like [Kondaurov and Nikitin, 1987; Takamizawa and Matsuda, 1990; Takamizawa, 1991]. Takamizawa and Matsuda 
[1990] and Takamizawa [1991] considered a local stress-free configuration by using a multiplicative decomposi- 
tion of deformation gradient, although this decomposition is implicit in their presentation. They realized that 
the local stress-free configurations are not unique but a corresponding metric is unique and defined a global 
stress-free configuration by equipping the original reference configuration with this metric giving it a Riemannian 
structure. Here we look at this metric and its rigorous connection with F = F e F g . 

It should be mentioned that similar ideas were used in plasticity and thcrmoclasticity before the growth 
mechanics applications. For the less familiar application in thermal stresses, the idea of decomposition of 
deformation gradient goes back to the works of Stojanovic and his coworkers [Stojanovic, et al., 1964; Stojanovic, 
1969]. See Vujosevic and Lubarda [2002]; Lubrada [2004] and Ozakin and Yavari [2009] for a detailed discussion. 
These researchers extended Kondo's [Kondo, 1955a,b, 1963, 1964] and Bilby's [Bilby, et al., 1955, 1957] idea of 
local elastic relaxation in the continuum theory of distributed defects to the case of thermal stresses. 13 

We have posed the following question in this paper: which space, as opposed to the Euclidean space, would 
be compatible with a relaxed state of the body? We claim that the answer to this question is, a Riemannian 
manifold whose metric is related to the nonuniform growth. This metric describes the relaxed state of the 
material with respect to which the strains in a given configuration should be measured. In this framework, 
the constitutive relations are given in terms of the material metric, the (Euclidean) spatial metric, and the 
deformation gradient F. 

Let us consider one of the above-mentioned imaginary relaxed pieces. Relaxation of this piece corresponds 
to a linear deformation (linear, because the piece is small) denoteed by F g . If this piece is deformed in some 
arbitrary way after the relaxation, one can calculate the induced stresses by using the tangent map of this 
deformation in the constitutive relations. In order to calculate the stresses induced for a given deformation of 
the full body, we focus our attention to one such particular piece. The deformation gradient of the full body at 



13 Note that the idea of local elastic relaxation was first proposed in the work of Eckart [1948]. 
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this piece F can be decomposed as F = F e F g , where, by definition, F e = FF~ 1 . Thus, as far as this piece is 
concerned, the deformation of the body consists of a relaxation, followed by a linear deformation given by F e . 
The stresses induced on this piece, for an arbitrary deformation of the body, can be calculated by substituting 
F e in the constitutive relations. One should note that F e and F g are not necessarily compatible. However, as 
long as we have a prescription for obtaining F e and F g directly for a given deformation map ip for the body and 
a growth distribution, we can calculate the stresses by the following procedure. Note also that if the material 
manifold is one dimensional the decomposition of deformation mapping into elastic and growth parts is always 
possible. This is implicitly assumed, for example, in [Scnan, et al., 2008]. 
For isotropic growth, one has the following expression for F g : 

(F 9 ) A b = S Sb' (3-1) 

Given this formula for F ff , we can calculate F e = FF" 1 for a given deformation, and use a constitutive relation 
that gives the stresses in terms of F e . At first glance these two approaches seem very different, however, they 
are related, as we demonstrate next. Wc should emphasize that the following discussion is not restricted to 
isotropic growth; given any F g our arguments can be repeated. 

The constitutive relations of the two approaches are formulated in terms of different quantities: G(X, t) and 
F on one side, and F e = FF~ 1 on the other. Let us start with our approach, namely, assume that a constitutive 
relation is given in terms of G(X, i) and F. This takes the form of a scalar free energy density function that 
depends on G(X, t), F, as well as on the spatial metric tensor g, and possibly X explicitly: 

* = *(X,e,G(X,t),F,goy)). (3.2) 

G, F, and g arc tensors, expressed in terms of specific bases for the material and the ambient spaces. A change 
of basis changes the components of these tensors, but iff does not change as it is a scalar. Let us consider a 
change of basis from the original coordinate basis E^ of the material space, with the following property 

((Ea,Eb))g = Gab , (3.3) 

to an orthonormal basis E^ that satisfies 

((E A ,E S )) G = S M . (3.4) 
A matrix F A B represents the transformation between the two bases: 

E A = F A B E B . (3.5) 

The orthonormality condition gives 

f A C ^B D GcD=S A6 . (3.6) 

Any F A c that satisfies this equation gives an orthonormal basis. Given such an F A C , we can also obtain an 
orthonormal basis for the dual space by using its inverse. Defining F c d as the transposed inverse of the matrix 
F A B , i.e., F A B F A c = 5 B and F A B F C b = we obtain the dual orthonormal basis {E" 4 } in terms of the original 
dual basis {E A } by 

E A = F A B -E B . (3.7) 
For isotropic growth, Gcd = e 2n( - x ^ScD = g(X, t) 2 5cD gives 14 

FA c = 6 c e -a^ t) = s c g -i {x ^ (39) 

as a solution to (3.6). Here, dh 3 is 1 for A = B, and 0, otherwise, i.e., 5^ = 8? = 5? = 1, etc. One should note 
that (3.6) has other solutions, as well, which we will comment on in the sequel. Now let us write the components 



14 This means that 

e 2n = g 2 or e n = g, (3.8) 
as g is always positive. If growth is anisotropic, having an expression for Gqd au these arguments can be repeated. 
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of the total deformation gradient F in the orthonormal basis {E^}. The components are transformed by using 
F as: 

pa A = ^A Bpa B ■ (3.10) 

Using (3.9), (3.8), and (3.1), we can clearly see that the components F a s arc given precisely by those of F e , 
the "elastic part" of the deformation gradient in F = F e F g approach: 

pa A = ^A Bpa B = Sfe- n ^F a B = (g(X, <))-Mf F a B = {F~ 1 )a B F a B = (F e ) a A . (3.11) 

Thus, F e is the original deformation gradient, written in terms of an orthonormal basis in the material space. 15 
We have also shown that there is no need for a mysterious "intermediate configuration" as the target space of 
F g , the latter simply gives an orthonormal frame in the material manifold, and as such, can be treated as a 
linear map from the tangent space of the material manifold to itself. 

Although a coordinate basis {E^ = d/dX A } is not necessarily orthonormal, one can always obtain an 
orthonormal basis by applying a pointwise change of basis ^ A B ■ Moreover, giving an orthonormal basis in 
this way is equivalent to giving a metric tensor at each point; the inner product of any two vectors can be 
calculated by using their components in the orthonormal basis. We have seen above that in the context of 
growth mechanics, this means that a change in the material metric due to a growth distribution can be given 
in terms of the "growth deformation gradient" of the local relaxation approach. Given an orthonormal basis 
{E^}, it is possible to obtain another one, {E^}, by using an orthogonal transformation A^ B : 

E' A = A/ E b , (3.12) 

where A^ B satisfies K. ^ K. g D 8 ^ ^ = 5 Ab - Let the relation between the original coordinate basis {E^} and the 
new orthonormal basis be given by the matrix F'^ 5 as follows 

E^ = F'/Eb. (3.13) 

The relation between F and F' is given as 



?' A B =k A °F e B . (3.14) 
Going in the opposite direction, one can see that F and F' represent the same material metric G, if and only if 



they are related through (3.14) for some orthogonal matrix A A B . This means that there is an 5*0(3) ambiguity 



in the choice of F, and hence, in that of F g . 

Using an orthonormal basis for the material manifold, we rewrite the constitutive relation (3.2) as 

* = *(X,6,Gab = 6 AB ,F a B = (F e ) a B ,g ab ). (3.15) 

Hence, given a constitutive relation <3> Rlcm in our (Ricmannian) approach, one can obtain a constitutive relation 
^ LR in the "local relaxation" approach by simply going to an orthonormal basis by (3.5) and (3.6), and ignoring 
the constant terms Gab = Sab and g ab = 5 ab in the functional dependence. 

* LR (X, 9, (F e f B ) = * Ricm (X, 6, Gab = 5 AB ,F a B = (F e ) a B,9ab = S ab ) . (3.16) 

Going in the opposite direction is also possible; starting with a free energy function for the F = F e F g approach, 
one can derive an equivalent free energy in the geometric approach. 16 



15 Note that given F^ s , material metric can be recovered as 

Gab = F c aF d b$c£,- 

16 A simple example can make this clearer. Let us assume that free energy density in the classical approach is fitrCe- In 
components this reads 

* = KCe) A sS A& = H (F A A F%) (t B B F b B ) S ab S A * = p {F- A F b B 5 ab ) (f^F^^) = ^ A F b B 8 ab G AB . 
Thus, ^ = /itr G C. 
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Balance of Mass. In the F = F e F 9 approach mass balance reads S m = + potr~L g , where L s = F 9 F~1. 
Usually, it is assumed that growth is density preserving [Lubarda and Hoger, 2002]. We show that the term 
trL g is equivalent to | trc (^jr)) where by trc we emphasize the G-dependence of the trace operator. Note 
that 

tr (§) = d -^T GAB = it ( F ^ F ^-) {^n B S dh ) = 2f\f A A = 2trL, 

Incompressibility. In the F = F e F g approach incompressibility is equivalent to J e = 1. In the geometric 
theory incompressibility means J = 1. These arc equivalent as is shown below: 



1 = J = W dct F = 1 dct (F a d F d c ) = det F a 6 = J e . (3.17) 

^det(FV^ fl S M ) 

Absolutely parallelizable manifolds and their connection with growth mechanics. Whenever defor- 
mation is coupled with other phenomena, e.g. plasticity, growth/remodeling, thermal expansion/contraction, 
etc. all one can hope for is to locally decouple the elastic deformations from the inelastic deformations. Many 
related works start from a decomposition of deformation gradient F = F e F Q , where F e is the elastic deforma- 
tion gradient and F a is the remaining local deformation or anelastic deformation gradient. Given an (inelastic) 
growth deformation gradient, a vector in the tangent space of X G B, i.e. W 6 TyJ5 is mapped to another 
vector W = F a W. Traditionally, these vectors are assumed to lie in the tangent bundle of an "intermediate 
configuration." In the literature, intermediate configuration is not clearly defined and at first glance it seems 
to be more or less mysterious. These are closely related to parallelizable manifolds (or absolutely parallelizable 
(AP) manifolds) [Eisenhart, 1926, 1927; Youssef and Sid- Ahmed, 2007; Wanas, 2008]. In an n-dimensional AP- 
manifold M, one starts with a field of n linearly independent vectors {E^)} that span the tangent vector at 
each point. We denote the components of E(^) by Ef^x. The dual vectors, i.e. the corresponding basis vectors 

for the cotangent space are denoted by {E^)} with components |Ej^|. Note that 

E< A) E( B) =^ and E< A) Ef A) = S J T . (3.18) 
One can equip M with a connection Tj K such that the basis vectors {E^)} are covariantly constant, i.e. 17 

E{ A) | J = 0. (3.19) 

Note that 

^(A)\JK - ^(A)\KJ = LJK^A) + T L K J~&\a)\L- (3.20) 

Therefore, (3.19) implies that 

TZ'ljk = 0, (3.21) 
i.e., M is fiat with respect to the connection Tj K . Note that 

<9E 7 

E (A)|J = q X j + r jRT E (A)- (3.22) 

Thus 

E^^T + ri* = 0. (3.23) 

Hence 

Tjk--Vj -gxK-Vi -QXK- (3-24) 



17 Equivalcntly, the tangent bundle is a trivial bundle, so that the associated principal bundle of linear frames has a section on 
M. 
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This connection has been used by many authors, e.g. by Bilby, et al. [1955] and Kondo [1955a] for dislocations, 
by Epstein and Elzanowski [2007] for material inhomogeneities, and by Stojanovic, et al. [1964] for thermal 
stresses. This connection is curvature-free by construction but has a non-vanishing torsion. 

For a growing body, in the local charts {^^} and {t/ 7 } for the reference and intermediate configurations, 
we have 

dU 1 = (Fg) 1 A dX A . (3.25) 
{Fg) 1 a can be identified with E^, and hence 

^ I jK = {F g ) I A J - (3.26) 

Note that this (growth) connection is curvature free but has a non-vanishing torsion. In plasticity it is shown 
that torsion of this connection has a physical meaning; it can be identified with the dislocation density tensor. 
For a growing body such a quantity does not seem to have a physical interpretation and we prefer to work with 
a Riemannian material manifold whose curvature quantifies the tendency of the growth distribution in causing 
residual stresses. 

In summary, our geometric approach has a concrete connection with that of F = F e F 9 : in the geometric 
approach we use a Riemannian manifold with a time-dependent metric as the material manifold while F = F e F s 
implicitly uses the same metric but in an absolutely parallelizablc manifold that is not Riemannian. We believe 
that our approach is more straightforward as we do not introduce an unnecessary torsion in the material manifold 
but of course the Riemannian material manifold has a non-vanishing curvature tensor, in general. 



4 Linearized Theory of Growth Mechanics 

Geometric linearization of elasticity was first introduced by Marsdcn and Hughes [1983] and was further de- 
veloped by Yavari and Ozakin [2008]. See also Mazzucato and Rachele [2006]. In this section, we start with 
a body with a time-dependent material manifold and its motion in an ambient space, which is assumed to be 
Euclidean. Suppose a given body with a material metric G is in a static equilibrium configuration, Lp. The 
balance of linear momentum for this material body reads 18 

DivP + /9 B = 0. (4.1) 

Now suppose the body grows by a small amount represented by a small change in the material metric SG. Lp will 
no longer describe a static equilibrium configuration. Stress in this new equilibrium configuration tp' = ip + Sip 
will be P' = P + <5P. We are interested in calculating the change in the stress (or the configuration), for a given 
small amount of growth. 

The linearization procedure can be formulated rigorously if instead of thinking about two nearby config- 
urations and the differences between various quantities for these configurations, we describe the situation in 
terms of a one-parameter family of configurations around a reference motion, and calculate the derivatives of 
various quantities with respect to the parameter. Let G e (X) be a one-parameter family of material metrics, 
ip e be the corresponding equilibrium configurations, and P e be the corresponding stresses. Let e = describe 
the reference configuration. Now, for a fixed point X in the material manifold, y e (X) describes a curve in the 
spatial manifold, and its derivative at e = gives a vector U(X) at <p(X) [Yavari and Ozakin, 2008]: 

= ^1 . (4>2) 
de e=o 

Considering Sip m t^f-, we see that a more rigorous version of Sep is the vector field U. U is the geometric 
analogue of what is called displacement field in classical linear elasticity. 
First variation (or linearization) of deformation gradient is defined as 



£(F) := V o F e 

8e 



V a 



* V ax 



= VU. (4.3) 

e=0 



18 Growth is a "slow" process compared to elastic deformations and hence inertial effects can be ignored. Throughout this paper, 
time is teated as a parameter. 
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Or in components 



C(F) a A = U a lA = 7 r ¥J + l a bc F b A U c , (4.4) 



dU a 

where 7? are the connection coefficients of the Ricmannian manifold (<S, g). Note that for different values of 

d_ 



e the spatial leg of F e lies in different tangent spaces and this is why covariant derivative with respect to c 



should be used. The right Cauchy-Green strain tensor for the perturbed motion tp t ^ is defined as 

Cab (c) = F a A (e)F b B (e)g ab (e). (4.5) 
Note that C e lies in the same linear space for all e € I, and the first variation of C can be calculated as 

^-C AB (e) = W JL F a A (e)F b B (e)g ab (e) + F a A (e)V jtF b B (e)g ab (e). (4.6) 

ae a « 3e 

Therefore 

C(C) AB := ^ C AB {e) = F b B g ab U a \ A + F a A g ab U b \ B . (4.7) 
ae e=o 

Transpose of the deformation gradient has the following linearization [Yavari and Ozakin, 2008]: C (F T ) = 
(VU) T . Spatial and material strain tensors are defined, respectively, as [Marsdcn and Hughes, 1983] 

e=i(g-pt.G) and E=i(*$g-G). (4.8) 

Or in components 

e ab = \ (g ab - G AB F- A a F~ B b ) , E AB = ^(C AB -G AB ). (4.9) 
We now show that linearization of E is related to e = ^£ u g, where u = U o ip^ 1 . We know that 

C(C) AB = g ab F a A F c B u b \ c + g ab F b B F c A u a w = F a A F C B u a{c + F b B F C A u b]c = 2F a A F b B e ab , (4.10) 
where e ab = ^(u a \ b + u b \ a ) is the linearized strain. Therefore 

£(C) = 2<p* t e. (4.11) 

Thus 

e = <p u C(E). (4.12) 

In other words, linearized strain is the push-forward of the linearized Lagrangian strain. Obviously, if the ambient 
space is Euclidean and the coordinates are Cartesian the covariant derivatives reduce to partial derivatives and 
one recovers the classical definition of linear strain in terms of partial derivatives, i.e. 

1 (du a du b \ 

Note that when the linearized strain is zero the variation field is a Killing vector field for the spatial metric g. 
In other words, this shows that this definition of linearized strain is consistent when the variation field generates 
an isometry of the ambient space. 

For the one-parameter family of material metrics G £ , variation of the material metric is defined as 



<5G w e — 

de 



G £ . (4.14) 

e=0 



In the case of isotropic growth 



d_ 



*j e = -re bo = z— — 

3 ae ae 



G = /3G, (4.15) 

e=0 



where f3 = 28 fl. Now consider, in the absence of body forces, the equilibrium equations DivP = for the 
family of material metrics parametrized by e : Div e P e = 0. Linearization of equilibrium equations is defined 
as [Yavari and Ozakin, 2008]: 

(Div,P £ )=0. (4.16) 



de 



e=0 
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Once again, one should note that since the equilibrium configuration is different for each e, P e is based at 
different points in the ambient space for different values of e, and in order to calculate the derivative with 
respect to e, one in general needs to use the connection (parallel transport) in the ambient space. For the case 
of Euclidean ambient space that we are considering and a Cartesian coordinate system {a; a }, (4.16) is simplified 
and in components reads 

dP aA (e) , . 



Thus, the linearized balance of linear momentum can be written as 



d d 



dX A de 



=o 



P aA {e) + 



r 



e=0 



AB \ 



Note that 



paA _ gdc 



paB , -pA u 

F ~de 



e=0 



P aB (e) = 0. 



dF°» 



(4.17) 

(4.18) 
(4.19) 



where $ = ^(X, 9, F, G, g) is the material free energy density. In calculating dP dc — , we need to consider the 
changes in both F and G: 



dP aA (e) _ dP aA dF b B dP aA dG CD 



Let us define 



A aA 



QpaJ 



de dF b B de 

d 2 ^ 



dG 



CD 



de 



b — 



dF b B 9 dF b sdF c a ~ G CD * dGcndF* A : 

where the derivatives are to be evaluated at the reference motion e = 0. Noting that for the case of an Euclidean 
ambient space (sec (4.4)) 

dF a A 

£=0 



and 



faACD 



paA 



a 2 * 



(4.20) 



(4.21) 



de 



dU a 



dX A 



we obtain 



Using 



and 



we obtain 



d_ 

de 



=o 



P aA {e) = K aA b B U b , B + B aACD SG C D. 



1 



AD 



[dG 



BD 



dG 



CD 



BC 



c 

2 V dX c ' 8X B 
dG AB 



dG 



BC 



ox D 



ACnBD 



de 



-G G 



dGci 
de 



5Ta * = de 



T A AB {e)=-G CD 8GcDT A A B + 7 :G AD 



dSG 



BD 



dSG 



CD 



dSG 



BC 



8X C 



dX 1 



dX 1 



In the case of isotropic growth, this is reduced to 

d_ 
de 



AB 



(«)=n 



3 dp 



e=o — ■ ' 2dX B ' 
With these results, the linearized balance of linear momentum (4.16) becomes 

3 df3 



(4.22) 
(4.23) 
(4.24) 
(4.25) 

(4.26) 
(4.27) 



(A aA b B U\ B ) A + {B aACD SG CD ) 



paB = g_ 



> A 2dX B 

Assuming that A and IB are independent of X, the linearized equilibrium equations are simplified to read 



(4.28) 



A B 



d 2 U b 



,aACD 



dX A dX B 

If the initial configuration is stress-free, we have 

r) 2 TT b 

j^aA.B ° U 



Gcd 



d(3 3 dp 



dX A 2dX B 

op 



pa 



B 



0. 



-B aACD G CD 



dX A dX B °^ dX A ' 

Let us now simplify the above linearized equations for a specific class of elastic materials. 



(4.29) 



(4.30) 
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Saint- Venant-KirchhofF materials. Saint- Venant-Kirchhoff materials have a constitutive relation that is 
analogous to the linear isortropic materials, namely, the second Piola-Kirchhoff stress S is given in terms of the 
Lagrangian strain E = i(C — G) as [Marsden and Hughes, 1983] S = A(trE)G -1 + or in components 

S CD = XE AB G AB G CD + 2fiE CD = ^{C AB G AB - 3)G CD + ^{C AB G AC G BD - G CD ), (4.31) 

where A = A(X) and fi = /u(X) are two scalars characterizing the material properties. We can obtain the tensor 
B aCAB from g as follow8 

>*CAB _ 9 / ab \ _ dP« C dS CD 



Using 

^— = -G AM G BN (4.33) 

we obtain 

B aACD G CD = —2C mn F a b (\G AB G MN + 2 l iG AM G BN ) + (3A + 2 f i)F a B G AB . (4.34) 

The initial metric is Euclidean; in Cartesian coordinates, Gab = Sab- Since the ambient space is also Euclidean, 
we can choose a Cartesian coordinate system whose axes coincide with the initial location of the material points 
along the material Cartesian axis. This will give, F a A — S A , where a and A both range over 1, 2, 3. Hence 

B aACD G CD = -^±^t 5 aA . (4.35) 
Similarly, for an initially stress-free material manifold, we obtain 

A aA b B = F a MF c Ngbc \^qAM qBN + ^qABqMN + qAN qBM^ ^ ^ 

For the case of an initially Euclidean material manifold with Cartesian coordinates we have 

K<lA b B QXA^XB = (A + M ' ab + ^ UaM - (4 ' 37) 

Therefore, Eq. (4.30) reads 

3A + 2a 88 , 

(x + ii)u btab + fju a , u = — y^^i ( 4 - 38 ) 

where we have identified the indices a and A. In analogy with thermal stresses, BS a b can be thought of as an 
cigenstrain. See Goriely, ct al. [2008] for a review of the existing linearized growth models. 

Stress-free growth distributions in the linearized theory. In this paragraph we show that in dimension 
three if [3 is linear in {AT^}, i.e. if j3 = a • X for some constant vector a, then a stress-free body remains 
stress free after growth. This is very similar to what is already known in classical linear thermoelasticity: 
temperature distributions linear in Cartesian coordinates leave a stress-free body stress free [Boley and Weiner, 
1997; Ozakin and Yavari, 2009]. 

Let us consider a one-parameter family of material metrics G £ and assume that the initial material metric 
is Euclidean, i.e. G £= o = 5. The corresponding curvature tensor is 7?. c . We need to calculate the linearized 
curvature, i.e. 

at 



e=0 



(4.39) 



This will give the solution to stress- free growth distributions. Note that SG = ^| e G e corresponds to a 
linearized growth and is stress-free if and only if 57L vanishes. To calculate the curvature variation, we follow 
Hamilton [1982] and denote derivative with respect to e by '. Following the definition of curvature tensor we 
can write 



i?' — 1 ( ® 2 G'bd 9 2 G' bc d 2 G' AD d 2 G' AC \ 1 nPQ , , 
K ABCD -- y dxAdx c ~ dX A dX D - dX B dX c + WBdxD ) + 2 G {KabcpG qd + K ABPD G 



(4.40) 
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In the case of Ricci curvature 

R'ab = G CD TZ' ACBD + {G CD ) TZacbd = G CD TZ' ACBD - G CP G dq G' pq TZ A cbd- (4-41) 
Similarly, for scalar curvature we have 

R' = g AB R' AB + (G AB )' Rab = G AB R' AB - G AP G B ^G' PQ R AB . (4.42) 
If the initial material manifold is Euclidean, i.e. if TZacbd = and Rab = 0, we have 

1 / d 2 5G BD d 2 5G BC d 2 8G AD OHGac \ . . 

0Kabcd 2 \dX A dX c 8X A dX D 8X B dX c dX B dX D )' [ ' 

1 / dHGcp dHG BC BHGad dHG AB \ , CD , . 

AB 2 \dx A dx B dx A dx D dx B dx c dx c dx D ) ' [ ' 

§ d 2 SG BC ABgCD _ dHG AB S AB § CD (4 45) 

dx A dx D dX c dX D • [ ' 

In the case of isotropic growth we have SGab — P&ab- In dimension three, vanishing of Ricci curvature is 
equivalent to vanishing the curvature tensor. Thus, SRab = reduces to 



d 2 P d 2 f3 xCD 



dX A dX B dX c dX D 

This is equivalent to 



Sab = 0. (4.46) 



A 12 = /? 13 = #23 = 0, (4.47) 

20,11 + 0,22+0,33 = 0, (4.48) 

0,11 + 20.22 +033 = 0, (4.49) 

0,11 + 022 + 2033 = 0. (4.50) 

The three relations (4.47) imply that = f{X v ) + g(X 2 ) + h(X 3 ) for arbitrary functions /, g, and h. The next 
three relations (4.48)-(4.50) imply that 0, u = 0, 22 = 033 = and therefore /"(X 1 ) = g"{X 2 ) = ti'(X 3 ) = 0, 
and hence is linear in Cartesian coordinates of the initial material manifold. 
In dimension two, SR = reduces to 

dx^ SAB - °' < 4 - 51 ) 

This means that has to be a harmonic function to represent a stress-free growth distribution. Again, this is 
very similar to what we know from classical linear thermoelasticity [Bolcy and Weiner, 1997; Ozakin and Yavari, 
20091. 



5 Concluding Remarks 

In this paper, we presented a geometric theory of elastic solids with bulk growth. We assumed that material 
points are preserved but density and "shape" are time dependent. We modeled a body with bulk growth by a 
Riemannian material manifold with an evolving metric tensor. The time dependency of material metric is such 
that the growing body is always stress free in the material manifold. We showed that energy balance needs to 
be modified when material metric is time dependent. Covariance of energy balance then gives all the balance 
laws. We also showed that entropy production inequality has a non-standard form when material manifold 
has an evolving metric. We showed that a more general notion of covariance of energy balance that includes 
temperature rescalings, in addition to giving all the balance laws, gives the constitutive restrictions imposed 
by the Clausius-Duhem inequality. We then showed how the principle of maximum entropy production can be 
used to obtain thermodynamically-consistent evolution equations for the material metric. 

We showed how analytical solutions for the residual stress field can be obtained in three examples of growing 
bodies with radial symmetries. We showed that even if mass is conserved, i.e. when growth results in only 
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shape changes, still one may see residual stresses. In the case of isotropic growth, we studied stress-free growth 
distributions using the material curvature tensor in both two and three dimensions. 

A concrete connection was made between our geometric theory and the conventional decomposition of 
deformation gradient into elastic and growth parts. We showed that in a special coordinate basis F e is our F. 
The present geometric theory is more natural and does not introduce a mysterious intermediate configuration. 
We linearized the nonlinear theory about a reference motion. Assuming that both the ambient space and the 
initial material manifold are Euclidean, we showed that growth results in eigen strains very similar to those 
of classical linear thermoelasticity. We found those growth distributions that are stress free in the linearized 
framework in both dimensions two and three. 
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